#include "quantum.h"

Qntm_var::Qntm_var(const Qntm_var & op) {
    
  int i;
    
  phi = op.phi;
  density = op.density;
  num = op.num;
  Eo = op.Eo;
  offset = op.offset;
  leftedge = op.leftedge;
  eval = new Eigenvalue [num];
  wave = new Fwave [num];
  prob = new double [num];
  if(!eval || !wave || !prob) nrerror("QNTM_VA::COPY CONSTRUCTOR: Memory allocation failure");

  for(i=0;i<num;i++) {
    eval[i]=op.eval[i];
    wave[i]=op.wave[i];
    prob[i]=op.prob[i];
  }

}

Eigenvalue Qntm_var::ret_eval(int i) {
  if(i>=num || !eval) nrerror("Qntm_var::ret_eval: error");
  return eval[i];
}
Fwave Qntm_var::ret_wave(int i) {
  if(i>=num || !wave) nrerror("Qntm_var::ret_wave: error");
  return wave[i];
}
Eigenvalue Qntm_var::ret_eval(double E) {
  int i = 0;
  while(i<num && eval[i].energy!=E) i++;
  if(i==num) nrerror("Qntm_var::ret_eval: error");
  return eval[i];
}
double Qntm_var::ret_energy(int i) {
  if(i>=num || eval==0) nrerror("Qntm_var::ret_energy: error");
  return eval[i].energy;
}
Fwave Qntm_var::ret_wave(double E) {
    int i = 0;
    while(i<num && eval[i].energy!=E) i++;
    if(i==num) nrerror("Qntm_var::ret_eval: error");
    return wave[i];
}

Qntm_var Qntm_var::operator=(Qntm_var & op) {
    int i;
    
    phi      = op.phi;
    density  = op.density;
    num      = op.num;
    Eo       = op.Eo;
    offset   = op.offset;
    leftedge = op.leftedge;
    eval     = new Eigenvalue [num];
    wave     = new Fwave [num];
    prob     = new double [num];
    if(!eval || !wave || !prob) nrerror("QNTM_VAR OPERATOR =: Memory allocation failure");

    for(i=0;i<num;i++) {
      eval[i]=op.eval[i];
      wave[i]=op.wave[i];
      if(op.prob) prob[i]=op.prob[i];
    }
    return *this;
}


void Qntm_var::show_wave(double E) {
  int i=check_norm(E);
  if(i>=0) wave[i].show();
  else printf("Show_wave: no wave with energy %le\n",E);
}


/* Create 3 arrays containing the eigenvalues, the wavefunctions and the probabilities */
/* Sort this arrays in energy order.It needs the array top with the list of eigenvalues */
/* 'energy' is the reference energy value, 'x_coor' is the position in the device in unit of length */
/* 'x_left' contains the starting position of the quantum structure. */
/* 'mass' is the relative mass in the mesh point */
void Qntm_var::set_eigen(Eigenvalue * top, double energy, double x_coor,double x_left, double mass) {
  short  swap;
  int 	 j;
  double buffer[4];
  Complex cbuffer;
  Eigenvalue *tmp;

    Eo 	     = energy;
    leftedge = x_left*LNTHCONV;
    offset   = (x_coor*LNTHCONV)-leftedge;
    tmp      = top;
    num	     = 0;
    while(tmp!=NULL) {
      num++;
      tmp = tmp->ptr;
    }
    if(eval) delete [] eval;
    if(wave) delete [] wave;
    if(prob) delete [] prob;
    eval = 0;
    wave = 0;
    prob = 0;
    if(num>0) {
      eval = new Eigenvalue [num];
      wave = new Fwave [num];
      prob = new double [num];
      j	 = 0;
      tmp  = top;
      for(j=0;j<num;j++) {
	eval[j].energy = tmp->energy;
	eval[j].sigma  = tmp->sigma;
	eval[j].gamma  = tmp->gamma;
	wave[j].kv     = cplx_sqrt(mass*KKK*(eval[j].energy-Eo));
	wave[j].mass   = mass;
	prob[j] = 0.0;
	tmp = tmp->ptr;
      }
      do {
	swap = 0;
	for(j=0;j<num-1;j++) {
	  if(eval[j].energy>eval[j+1].energy) {
	    buffer[0] = eval[j].energy;
	    buffer[1] = eval[j].sigma;
	    buffer[2] = eval[j].gamma;
	    buffer[3] = wave[j].mass;
	    cbuffer   = wave[j].kv;
	    eval[j].energy = eval[j+1].energy;
	    eval[j].sigma  = eval[j+1].sigma;
	    eval[j].gamma  = eval[j+1].gamma;
	    wave[j].kv     = wave[j+1].kv;
	    wave[j].mass   = wave[j+1].mass;
	    eval[j+1].energy= buffer[0];
	    eval[j+1].sigma = buffer[1];
	    eval[j+1].gamma = buffer[2];
	    wave[j+1].kv    = cbuffer;
	    wave[j+1].mass  = buffer[3];
	    swap = 1;
	  }
	}
      } while(swap);
    }
}


/* Remove an eigenvalue */
int Qntm_var::drop(double E) {

  int j,i,ckr;
  Eigenvalue *tmp1;
  Fwave * tmp2;
  double * tmp3;

  ckr = 0;
  for(j=0;j<num;j++) if(eval[j].energy==E) ckr = 1;
  if(ckr==0) return 0;
  if(ckr==1 && num>1) {
    tmp1 = new Eigenvalue [num-1];
    tmp2 = new Fwave [num-1];
    tmp3 = new double [num-1];
    
    i = 0;
    for(j=0;j<num;j++) {
      if(eval[j].energy!=E) {
	tmp1[i] = eval[j];
	tmp2[i] = wave[j];
	tmp3[i] = prob[j];
	i++;
      }
    }

    if(eval) delete [] eval;
    if(wave) delete [] wave;
    if(prob) delete [] prob;
    eval = tmp1;
    wave = tmp2;
    prob = tmp3;
    num = num-1;
    return 1;
  }
  else {
    if(eval) delete [] eval;
    if(wave) delete [] wave;
    if(prob) delete [] prob;
    num = 0;
    return 1;
  }

}


int Qntm_var::set_wave(double E, Complex Ac, Complex Bc) {
    int i = check_norm(E);
    if(i>=0 && i<num) {
      wave[i].A = Ac;
      wave[i].B = Bc;
      return 1;
    }
    else return 0;
  }

/* ATTENTION: The field lgth has to be an offset-length expressed in cm */ 
double Qntm_var::calc_prob(double E, double lgth) {
    int i = check_norm(E);
    double x,sum = 0,h;
    if(i==-1) return 0;
    lgth = lgth*LNTHCONV;
    x =	offset;
    if(lgth==0) return cplx_mod_sqr(wave[i].eval(x));
    if(lgth<0) return cplx_mod_sqr(wave[i].eval(x-lgth));
    h = lgth/50.0;
    while((x-offset)<lgth){
      sum += cplx_mod_sqr(wave[i].eval(x));
      x += h;
    }
    return sum*h;
}



void Qntm_var::set_prob(double E, double lgth) {
    int i = check_norm(E);
    double x,sum = 0,h;
    if(i!=-1) {
      x = offset;
      h = lgth*LNTHCONV/50.0;
      while((x-offset)<lgth*LNTHCONV) {
	sum += cplx_mod_sqr(wave[i].eval(x));
	x   += h;
      }
      prob[i] = sum*h;
      if(prob[i]>1.0 || prob[i]<1e-20) prob[i] = 0;
    }
}


void Qntm_var::set_prob_e(double E, double lgth, FILE * fp) {
    int j,i = check_norm(E);
    double x,sum = 0,val,h;
    if(i!=-1) {
      x = offset;
      h = (LNTHCONV*lgth)/50.0;
      for(j=0;j<50;j++) {
	val = cplx_mod_sqr(wave[i].eval(x));
	fprintf(fp,"%le  %le\n",(x+leftedge)/LNTHCONV,val+E);
	sum += val;
	x   += h;
      }
      prob[i] = sum*h;
      if(prob[i]>1.0 || prob[i]<1e-20) prob[i] = 0;
    }
}

void Qntm_var::set_prob_h(double E, double lgth, FILE * fp) {
    int j,i = check_norm(E);
    double x,sum = 0,val,h;
    if(i!=-1) {
      x = offset;
      h = (LNTHCONV*lgth)/50.0;
      for(j=0;j<50;j++) {
	val = cplx_mod_sqr(wave[i].eval(x));
	fprintf(fp,"%le  %le\n",(x+leftedge)/LNTHCONV,-val-E);
	sum += val;
	x   += h;
      }
      prob[i] = sum*h;
      if(prob[i]>1.0 || prob[i]<1e-20) prob[i] = 0;

    }
}


  /* Check if the input energy is an exisiting eigenvalue energy */
  /* If that is the case the function returns the position index, else it returns -1 */
int Qntm_var::check_norm(double E) {
    int i=0;
    while(i<num && fabs(E-eval[i].energy)>1e-4) i++;
    if( i<num ) return i;
    else return -1;
  }


  /* This routine set the new wave function coefficients Ac and Bc between two energy steps */
  /* It multiplies the old ones with the transfer matrix to get the new ones. */
void Qntm_var::propagation(double E, Complex & Ac, Complex & Bc, Qntm_var * old) {

    int i = check_norm(E);
    int j = old->check_norm(E);
    Complex A_old,B_old,ii(0,1),tmp1,tmp2;

    if(i>=0 && j>=0) {
      Tmatrix T(old->wave[j].kv,wave[i].kv,old->wave[j].mass,wave[i].mass,ret_offset());
      A_old = Ac;
      B_old = Bc;
      Ac = T.prod1(A_old,B_old);
      Bc = T.prod2(A_old,B_old);

      /* checking if there is effective continuity
      coeff1 = cplx_mod((old->wave[j].mass*MASS/PLANKUT)/old->wave[j].kv);
      coeff2 = cplx_mod((wave[i].mass*MASS/PLANKUT)/wave[i].kv);
      tmp1 = A_old*coeff1*cplx_exp(old->wave[j].kv*ret_offset()*ii)
	+B_old*coeff1*cplx_exp(old->wave[j].kv*ii*(-ret_offset()));
      tmp2 = Ac*coeff2*cplx_exp(wave[i].kv*ret_offset()*ii)
	+Bc*coeff2*cplx_exp(wave[i].kv*ii*(-ret_offset()));
      if(fabs(tmp1.rl-tmp2.rl)>1e-3 ||fabs(tmp1.im-tmp2.im)>1e-3 ) {
	printf("Propagation error: %lf-%lf ... %lf-%lf\n",tmp1.rl,tmp1.im,tmp2.rl,tmp2.im);
	exit(1); 
	}*/

    }
    else {
      printf("Propagtion failure i=%d j=%d\n",i,j);
      show();
      old->show();
      Ac.rl = 0;
      Ac.im = 0;
      Bc.rl = 0;
      Bc.im = 0;
    }

  }

/* ------------------------------------------ */
/* ------------------------------------------ */
/* it returns 0 if the eigenvalue is already present */

int Qntm_var::add_eigenvalue(Eigenvalue eval_inp, Fwave wave_inp) {

  int i;
  double * tmp3;
  Eigenvalue * tmp1;
  Fwave * tmp2;

  for(i=0;i<num;i++) if(eval[i].energy==eval_inp.energy) return 0;

  num++;
  tmp1 = new Eigenvalue[num];
  tmp2 = new Fwave[num];
  tmp3 = new double[num];

  i = 0;
  while(i<num-1 && eval_inp.energy>eval[i].energy) {
    tmp1[i]=eval[i];
    tmp2[i]=wave[i];
    tmp3[i]=prob[i];
    i++;
  }
  tmp1[i]=eval_inp;
  tmp2[i]=wave_inp;
  tmp3[i]=0;
  i++;
  while(i<num) {
    tmp1[i]=eval[i-1];
    tmp2[i]=wave[i-1];
    tmp3[i]=prob[i-1];
    i++;
  }

  if(eval) delete [] eval;
  if(wave) delete [] wave;
  if(prob) delete [] prob;

  eval = tmp1;
  wave = tmp2;
  prob = tmp3;
  
  return 1;

}
/* ------------------------------------------ */
/* ------------------------------------------ */
/* it returns 0 if the eigenvalue is already present */

int Qntm_var::add_eigenvalue(double E) {

  int i;
  double * tmp3;
  Eigenvalue * tmp1;
  Fwave * tmp2,tmp4;

  if(num==0) nrerror("ADD_EIGENVALUE: No other eigenvalues!");
  for(i=0;i<num;i++) if(eval[i].energy==E) return 0;

  num++;
  tmp1 = new Eigenvalue[num];
  tmp2 = new Fwave[num];
  tmp3 = new double[num];

  tmp4 = wave[0];

  i = 0;
  while(i<num-1 && E>eval[i].energy) {
    tmp1[i]=eval[i];
    tmp2[i]=wave[i];
    tmp3[i]=prob[i];
    i++;
  }
  tmp1[i].energy=E;
  tmp1[i].sigma=0.0;
  tmp1[i].gamma=0.0;
  tmp1[i].ptr = 0;
  tmp2[i]=tmp4;
  tmp2[i].kv=cplx_sqrt(tmp2[i].mass*KKK*(E-Eo));
  tmp3[i]=0;
  i++;
  while(i<num) {
    tmp1[i]=eval[i-1];
    tmp2[i]=wave[i-1];
    tmp3[i]=prob[i-1];
    i++;
  }

  if(eval) delete [] eval;
  if(wave) delete [] wave;
  if(prob) delete [] prob;

  eval = tmp1;
  wave = tmp2;
  prob = tmp3;
  
  return 1;

}

/* ------------------------------------------ */
/* ------------------------------------------ */
/* it returns 0 if the eigenvalue is already present */

int Qntm_var::add_eigenvalue(double E, double mass_inp) {

  int i;
  double * tmp3;
  Eigenvalue * tmp1;
  Fwave * tmp2;

  for(i=0;i<num;i++) if(eval[i].energy==E) return 0;

  num++;

  tmp1 = new Eigenvalue[num];
  tmp2 = new Fwave[num];
  tmp3 = new double[num];
  


  i = 0;
  while(i<num-1 && E>eval[i].energy) {
    tmp1[i]=eval[i];
    tmp2[i]=wave[i];
    tmp3[i]=prob[i];
    i++;
  }
  tmp1[i].energy=E;
  tmp1[i].sigma=0.0;
  tmp1[i].gamma=0.0;
  tmp1[i].ptr = 0;
  tmp2[i].mass = mass_inp;
  tmp2[i].kv=cplx_sqrt(tmp2[i].mass*KKK*(E-Eo));
  tmp3[i]=0;
  i++;
  while(i<num) {
    tmp1[i]=eval[i-1];
    tmp2[i]=wave[i-1];
    tmp3[i]=prob[i-1];
    i++;
  }

  if(eval && num>1) delete [] eval;
  if(wave && num>1) delete [] wave;
  if(prob && num>1) delete [] prob;

  eval = tmp1;
  wave = tmp2;
  prob = tmp3;

  return 1;

}

/* ------------------------------------------ */
/* ------------------------------------------ */

void Qntm_var::modify_wave(Eigenvalue eval_inp, Fwave wave_inp) {

  int i;

  i = 0;
  while(i<num && eval[i].energy!=eval_inp.energy) i++;
  if(i>=num) printf("MODIFY WAVE: eigenvalue not found!\n");
  else wave[i] = wave_inp;

}

/* ------------------------------------------ */
/* ------------------------------------------ */

void Qntm_var::show() {

  int i;
  printf("Num=%d  Eo=%le  leftedge=%le  offset=%le\n",num,Eo,leftedge,offset);
  printf("Phi=%le  density=%le\n",phi,density);
  for(i=0;i<num;i++) {
    printf("%d Eigenvalue:\n",i);
    eval[i].show();
    printf("%d Wave:\n",i);
    wave[i].show();
    printf("%d Prob:\n",i);
    printf("%le\n",prob[i]);
  }
  printf("\n");
}










